Two-body correlations in A^-body boson systems 
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We formulate a method to study two-body correlations in a system of A*' identical bosons inter- 
acting via central two-body potentials. We use the adiabatic hyperspherical approach and assume a 
Faddeev-like decomposition of the wave function. For a fixed hyperradius we derive variationally an 
optimal integro-differential equation for hyperangular eigenvalue and wave function. This equation 
reduces substantially by assuming the interaction range much smaller than the size of the A'^-body 
system. At most one-dimensional integrals then remain. We view a Bose-Einstein condensate pic- 
torially as a structure in the landscape of the potential given as a function of the one-dimensional 
hyperradial coordinate. The quantum states of the condensate can be located in one of the two 
potential minima. We derive and discuss properties of the solutions and illustrate with numerical 
results. The correlations lower the interaction energy substantially. The new multi-body Efimov 
states are solutions independent of details of the two-body potential. We compare with mean-field 
results and available experimental data. 

PACS numbers: 31.15.Ja, 05.30.Jp, 21.65.-|-f 



I. INTRODUCTION 

The average properties of an A^-body system are often 
investigated in the mean-field approximation ^ 
where the wave function is a product of one-particle am- 
phtudes. This excludes a priori effects of particle corre- 
lations, which often are responsible for decisive features 
or phenomena. A prominent example is the use of the 
bare nucleon-nucleon interaction in computations of nu- 
clear ground state structure. The Hartree-Fock results 
are catastrophic, either producing unbound systems or 
collapsed point-like structures with infinite binding en- 
ergy Q. The correlations must be dealt with here. One 
way is to use effective interactions and maintain the same 
Hilbert space of independent particles jsj. One can get 
a long way by this procedure, but the particles are still 
not correlated in the wave functions, although reasonable 
sizes and binding energies are found j|, |[ ^. Any effect 
of correlations can therefore not be tested experimentally 
except as deviations from the mean- field results. 

Another example, where correlations are essential, is 
the decay of Bose-Einstein condensates j6[ , which can- 
not be referred to as the lowest energy solution. Indeed, 
the A'^-body system has many lower lying states as imme- 
diately realized from the fact that two and more atoms 
form bound states. The energy is then already lower than 
the energy of the condensate, which therefore eventually 
decays into these lower lying structures. The decay can 
be either three-body recombination possibly enhanced by 
the presence of other particles 0, ||, or macroscopic 
collapse into a structure of very small spatial dimension, 
which subsequently then recombine into more favorable 
bound states ||, 

The direct recombination process is increasingly prob- 
able with larger scattering length |6[ and also the 
macroscopic collapse must increase strongly with scat- 
tering length . The limit of infinite scattering length 
is traditionally considered as difficult to solve as exem- 



plified by the three-body system where the delicate Efi- 
mov states could occur The structure of the 
Af-body system in this limit is at least as difficult as 
the three-body problem ||l5| , but the growing interest 
]16|, [7[ m and the difficulties in this area demand 
new approaches. 

The macroscopic collapse is conceptually very similar 
to the nuclear fission process where the (liquid) nucleus 
in a collective process is divided into two or more pieces. 
One difference is that in fission the fragments move away 
from each other, whereas (gaseous) Bose-Einstein con- 
densates first collapse into a more dense state and af- 
ter recombination into smaller subsystems the fragments 
move apart | pT| . These collapse mechanisms can to a 
large extent be described in a mean- field picture ^ , 
but especially for condensates the correlations could be 
very important first to establish the collective coordinate 
and the corresponding potential energy and second by 
influencing the collapse process itself . 

The essential ingredients in a description of correla- 
tions in A^-body systems are the techniques used to solve 
few-body problems Isolated two-body systems 

are easily solved and the key is very likely in handling 
of the three-body problem. This expectation arises since 
two particles plus all remaining particles effectively is a 
three-body problem, and this type of two-body correla- 
tions beyond the mean-field inside the A^-body system is 
probably dominating. This is also the philosophy in the 
Faddeev ||2^ and Yakubovsky equations, where the 
two-body amplitudes eventually are the basic quantities. 
The Yakubovsky reduction is rigorous, but very cumber- 
some in its full glory. However, such formulations may 
provide inspiration to practical approximations. 

Unfortunately, it is impossible to include correlations 
directly on top of mean-field calculations with attractive 
zero range interactions [^5[ The Gross-Pitaevskii 

equation without a confining external trap can only be 
used for repulsion in a product wave function. The low- 
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est state for an attraction would correspond to a diver- 
gent collapsed non-physical solution, but an additional 
external field is able to hold metastable solutions at 
larger distances. Skyrme Hartree-Fock is more sophisti- 
cated and an attraction leading to a bound system does 
avoid collapsed solution s |5[| . Including correlations trig- 
gers the Thomas effect ]27| |, where collapsed three-b ody 
states would appear inside the many-body system p6| . 
Renormalization can be invoked to cure these divergences 
and maintain the simplicity of the zero-range interaction 
p8| , |2^ , but correlations are still not included. Obviously 
a finite range interaction would also prevent the disaster 
when correlations are allowed. 

To study correlations the form of the wave function 
must be flexible enough to include the corresponding de- 
grees of freedom. As we learned from the mean-field ap- 
proximation sizes and binding energies may be rather 
accurate even with an imprecise wave function pl| , pO| . 
Thus for a specific purpose it seems possible to design a 
relatively simple wave function, which for other purposes 
may be a rather poor approximation, but precisely ac- 
counting for the desired degrees of freedom. A promising 
form of a correlated wave function suggested for nucleons 
pl| was recently extended to more general systems [|2| . 
The simplest structure is clearly found in Bose-Einstcin 
condensates. Here was recently introduced a formula- 
tion with generalized hyperspherical coordinates and an 
adiabatic expansion with the hyperradius as the adia- 
batic coordinate pO| . Only the crude approximations 
of a zero-range interaction and the lowest (constant and 
thus non-correlated) hyperspherical angular wave func- 
tion were used in |£0| . 

In this paper we shall use the simplifications of the 
completely symmetric wave function. These systems re- 
veal features only explainable by correlations and we can 
then already here obtain results of practical interest. The 
purpose of this paper is to formulate the details of a re- 



cently suggested |21| method to include two-body cor- 
relations in an iV-body system of identical bosons. We 
shall describe the qualitative features of the key quanti- 
ties, give numerical examples, compare with established 
methods and available experimental data and provide 
predictions going beyond the present knowledge. The 
model is apparently rich with many unexplored possibil- 
ities perhaps especially by application to systems with 
large scattering length, to the dynamical evolution, and 
to other systems of different symmetries. 



II. THE iV-BODY PROBLEM AND THE 
HYPERSPHERICAL FORMULATION 

A formulation accounting for particle correlations in an 
A^-body system must by definition go beyond the mean- 
field approximation. The Faddeev-Yakubovsky equa- 
tions could be an appropriate starting point. However, 
first the form of the Hamiltonian describing the system 
must be decided. Then a convenient set of coordinates 



must be chosen in harmony with the formulation and 
the anticipated approximations. Derivations of suitable 
equations of motion are then possible and their prop- 
erties can be investigated. This section describes how 
we choose two-body interactions, hyperspherical coordi- 
nates, a Faddeev-like decomposition of the wave function, 
use an adiabatic expansion, and assume at first only s- 
waves. 



A. The hyperspherical coordinates 

The system of N identical interacting bosons of mass 
m may be described by N coordinate vectors r; and mo- 
menta Pi, labeling the particles by the index i = 1, . . . , N . 
Here a more suitable choice of coordinates is the center 



of mass coordinates R — X^iLi ''^if-^j the A^ 
Jacobi vectors fjk with k = 1,2, . . . , N ~ 1 |32 
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and their associated momenta. These Jacobi coordinates 
are illustrated for up to six particles in fig. We use 
the notation ijk = \ffk\, so r/N^i is proportional to the 
distance between particles 1 and 2, 7]n-2 is proportional 
to the distance between particle 3 and the center of mass 
of 1 and 2, jyAr-a is proportional to the distance between 
particle 4 and the center of mass of the first three parti- 
cles, etc. 

Hyperspherical coordinates are now defined in relation 
to the Jacobi vectors. One length, the hyperradius p, is 
defined by 



2 _ \ ^ 2 
Pi =2^Vk 
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1 



N 

2 _ 2 \ ^ 

P =Pn-i = j^2^^' 

i<j 
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where r,- 



\f^-rj\. The A^-2 hyperangles Qifc G [0, 7r/2] 
for k — 2,3, . . . , N — 1 relate the length of the Jacobi 
vectors to the hyperradius via the definition 



sm ak 



Pk 



(3) 



Since pi — rji the variable ai — 7r/2 is superfluous, but 
is for convenience often included in the notation. Re- 
maining are the 2(A — 1) angles fi^,'^' = {9k, fk) for 
k — 1,2, . . . , N — 1 defining the directions of the A^ — 1 
rffe-vectors, i.e. 9k G [0,7r] and ipk € [0, 27r]. All an- 
gles are collectively denoted by = {ak,9k,'fk}, k — 
1,2, . . . ,N - 1. In total n and p amount to 3{N - 1) 
degrees of freedom and the center of mass coordinates R 
amount to three. These coordinates are connected by 



N N , 

Vr' = -Vr2 +-( 

, — 1 ^ 



N 



NR^ . (4) 
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The total volume element is 11^=1 ^^^i — 
N^/'^d'^RY\^~^ d'^ffk, where the part depending on 
relative coordinates is Y\k=i '^^^k- In hyperspherical 
coordinates this relative part becomes 



where Pr = J^i Pi is the total momentum and M = Nm 
is the total mass of the system, i.e. 



H = H, 



total -^cm 



^2 p2 

^2m 2M 



N-l 



3iV-4 



dn 



k=l 

dClk 



dni^'> dni'''> dn 



k-l 



di}'^^ = dak sin^ ak cos 



3fe-4 



ak 



(5) 

(6) 
(7) 



where dfl 



(k) 



dOkSinOkdipk is the familiar angular vol- 
ume element in spherical coordinates. The recursion 
stops at dfli — dfl^\ Since the angle un-i is re- 
lated directly to a two-body distance, ri2, by sinajv_i = 
riN-i/ Pn-1 = 1^12/ (V^p), the volume element in eq. (|^) 
related to this angle is especially important: dJli^"^'' = 

rfcKW-l sin^ CtAT-i COS"^^"*" Q!Ar_i. 

The angular volume integrals can be computed [p4| , 
i.e. / dn'-a'' = 0Fr(3(fc-l)/2)/(4r(3A:/2)) and / dd,^^ 
47r, where r(a:) is the gamma function. 

Using eq. (|^) an angular matrix element of an operator 
O for fixed p is then 

(*|d|$)f2- J dnN-i-i'*{p,n)d^p,n) , (8) 

which in general is a function of p. 



B. Hamiltonian 

We consider N identical particles of mass m interact- 
ing through short range two-body potentials. For com- 
pleteness we will throughout the paper include an ex- 
ternal trap confining the particles to a limited region of 
space through a simple three-dimensional harmonic oscil- 
lator potential muj^rf/2. The trap is relevant for Bose- 
Einstein condensation but can easily be omitted. With a 
two-body central potential Vij = V{rij) the total Hamil- 
tonian is given by 



total 



^ / «2 1 

^\2m 2 



9 ? 

muj rf 



N 



(9) 



i<3 



which, using eq. (Q), is separable into a part only in- 
volving the center of mass coordinates and a part only 
involving relative coordinates. We can see this by sub- 
tracting 



p2 

2M 



1 



(10) 



N N 

J2lmu;\f-^Mu;'R'+Y^V.,. (11) 



i<j 



Using eq. (^ and denoting the intrinsic kinetic energy 
operator by T = Pi I (2m) — (2M) we can write 



1 ^ 
H = f+ -muj^p^ + Vi 



(12) 



In hyperspherical coordinates T can be rewritten as p2, 



T 



2m 



9 „37V^4 9 



p2 



(13) 



with the dimensionless angular kinetic energy operator 
A^_j^ recursively defined by 



cos^ ak 



sin^ ak 



nf = - 



dc 



3fc — 6 — {3k — 2) cos 2ak d 
sin 2ak dak 



(14) 



(15) 



where hlk is the angular momentum operator associated 
with ffk- The recursion stops at — If. The angular 
kinetic energy operator is thus a sum of derivatives with 
respect to the various hyperspherical angles. Convenient 
transformations to get rid of first derivatives in eqs. ( |l^ ) 
and (|l|) are 



h^^'-P dp^ dp 

~(3W-4)/2 (d^ _ (3jV-4)(3iV-6) ^ ^(3Ar-4)/2 



9p2 



4p2 



n2 = sin ^ ak cos' 



-(3fe-4)/2 



ak 



92 



9fc- 10 



(3/c-4)(3fc-6) 



tan2 ttfc^ sin at cos^^'' **'/2 ^ (27) 



The Hamiltonian H can now be rewritten as 



1 



hn = A 



2 2 
muj p 



N 



2mp' 



i<j 



(18) 
(19) 



where Vij = 2mp^Vij /h^ is a dimensionless potential, Tp 
is the radial kinetic energy operator, and /iq is the dimen- 
sionless angular Hamiltonian. The intrinsic Hamiltonian, 
H, thus contains a part depending on p and a part, /in, 
depending on both p (parametrically) and fi. 
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C. Equations of motion 

Since the total Hamiltonian is given as -fftotai = ^cm + 
H, the total wave function for the TV-particle system can 
without loss of generality be written as a product of a 
function, T, depending only on R and a function, ^P, 
depending on p and the 3A^— 4 angular degrees of freedom 
collected in Vt: T{R)^{p,Q,). The center of mass motion 
for the total mass M = Nm is given by 



i/cmT(i?) = ^cmT(i?) 



(20) 



and as seen from eq. ( ]lO| ) the energy spectrum is that of 
a harmonic oscillator, i.e. Ecm.n = hLu{n + 3/2), where n 
is a non-negative integer. 

The relative wave function ^'(p, il), obeying the sta- 
tionary Schrodinger equation, 

H^ip.Q) ^ E'<S/ip,n) , (21) 

is for each value of the hyperradius p expanded as 

oo 

^{p, n) = p-(3A^-4)/2 ^ /„(p)$„(p, n) , (22) 

n=0 

I 



where the factor p~(3A'-4)/2 jg included to eliminate first 
derivatives in p, see eq. (|l6|). Here the hyperradial wave 
functions, fnip), are the expansion coefficients for fixed 
p on the complete set of solutions ^n{p, ^) obtained by 
solving the angular eigenvalue equation: 



(/in-A„)$„(p,r!) = 0, 



(23) 



where A„ is the angular eigenvalue, which depends on 
p. We will usually apply the normalization ($„|$m)o = 



In complete analogy to the technique employed for 
= 3 [|5) we insert eq. ( p2| ) in eq. (|2l]), use eqs. ( |l^ ) 
and (|2^), and finally project the resulting equation onto 
the angular eigenfunctions <&„(p, fl). We then arrive at a 
set of coupled radial equations 



2mE A„(p) (3jV-4)(3iV-6) p' \ ^ / (i) d (2) \ 



(P) , (24) 



where the trap length is bt 
terms Q^H (QiM 



h/ (muj) 
0) are defined as 



{'^n{p,n)\ 








|$„(P,^^))^ 



and the coupling 



(25) 



The angular eigenvalues A„ enter these coupled equations 
as a radial potential. The total diagonal effective radial 
potential, Un (p) , entering on the left hand side of eq. (|2^ ) 
is: 



2mf/„ _ A„ (37V - 4) (37V - 6) 



^4 ^nn 



(26) 



This includes a p^-term due to the external harmonic 
field, a p~^ centrifugal barrier-term due to the transfor- 
mation of the radial kinetic energy operator, the angular 

(2) 

potential A„, and the diagonal coupling term Qnn- 

If the non-diagonal coupling terms are neglected, 
i.e. the right hand side of eq. ( p4| ) vanishes, the equa- 
tions simplify significantly to 



2m dp^ 



E„ 



■n,q jjn,q 



fnAp) = . (27) 



We have here included another index q to indicate a spec- 
trum of energies and radial wave functions obtained for 
each angular solution n. 



D. The angular eigenvalue equation 



The eigenvalue A from eq. (|2^) is the key quantity car- 
rying essentially all information about the two-body in- 
teractions and therefore about possible correlations as 
well. The technique and approximations used to find A 
are then especially important. 



1. Decomposition of the angular wave function 

The angular eigenvalue is obtained by solving eq. (|22 
For each p we first assume a decomposition of the wave 
function $ (omitting the index n) in additive components 



I.e. 



N 



<I>(p,17) = ^$,,(p,l]), 



(28) 



i<3 



where each term $ij is a function of p and all angular 
coordinates f2. This decomposition is in principle exact, 
since each term in itself is sufficient when all 51 degrees of 
freedom are allowed. At first this ansatz seems clumsy by 
introducing an overcomplete basis. However, the indices 
i and j indicate special emphasis on the particle pair i—j. 
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The component $y is expected to carry the information 
associated with two-body correlations of this particular 
pair. This has no significance before it is exploited in 
numerical techniques or subsequent approximations. 

Rewriting the wave function obeying the Schrodinger 
equation as a sum of terms has been very successful in 
three-body computations. The advantage is that the cor- 
rect boundary conditions are simpler to incorporate as 
expressed in the original formulation by Faddeev in- 
tended for scattering. Still mathematically nothing is 
gained or lost in this Faddeev-type of decomposition. For 
very weakly bound and spatially very extended three- 
body systems, s-waves in each of the Faddeev compo- 
nents are sufficient to describe the system. This is ex- 
ceedingly pronounced for large scattering lengths where 
the delicate Efimov states appear [|l5[ ^ . 

The present A^-body problem is of course in general 
more complicated. However, for dilute condensates es- 
sential similarities remain, i.e. the relative motion of two 
particles on average far from each other is most likely 
dominated by s-wave contributions. Each particle can- 
not detect any directional preference arising from higher 
partial waves. Only the monopole prevails. Implement- 
ing these ideas in the present context imply that each 
amplitude $ij for a fixed p only should depend on the 
distance between the two particles. For that purpose 
we define a two-index parameter ay by 



sm ai 



(29) 



which is distinctively different from the a^'s of eq. (|^). 
Thus we assume 



(30) 



The boson symmetry implies that all the functions 
are equal and that we should not distinguish, so we there- 
fore omit the indices. We then arrive at the angular wave 
function 



N 
i<j 



N 

i<j 



(31) 



where we used (j>ij (p, ) = 4>{p, Uij ) = (l){aij ) with omis- 
sion of the coordinate p in the last notation. The wave 
function in eq. (31) is symmetric with respect to inter- 
change of two particles, i <-> j , since = aji and since 
terms like 4>{aik) + (t>{ajk) always appear symmetrically. 

This ansatz of only s-waves dramatically simplifies the 
angular wave function. The original overcomplete Hilbert 
space is now reduced, so not every angular wave function 
can be expressed in this remaining basis. Thus rigorously 
the reduction resulted in an incomplete basis, but the 
degrees of freedom remaining in eq. (|3^) are expected to 
be precisely those needed to describe the main features 
of the condensate. The approximations are tailored to 
the problem under investigation. 



S. Faddeev-like equations 

Inserting the ansatz of eq. (|3^) along with eq. ([l^) in 
eq. (E3h yields the angular equation 



N 



N 



(32) 



k<l 



with 
to 



(puj). Rearrangement of summations leads 



N 

E 



A ) (t)ij + V. 



N 

y E' 

k<l 



(33) 



For three particles the Faddeev equations are obtained 
by assuming that each term in the square brackets sep- 
arately is zero [35 . The same assumption for the iV- 



particle system results in the N{N — l)/2 Faddeev-like 
equations 



N 



E '^'^^ = ' 



(34) 



k<l 



which are actually identical due to symmetry. We shall 
not in the present paper rely on the validity of this as- 
sumption, but only use it to illustrate the procedure in 
the general discussion. 

Choosing i = \ and j = 2, with the ansatz for the 
wave function, eq. (|3l|), the kinetic energy operator A^_]^ 

in eq. ( p^ reduce to because A^_20i2 = and 

^Ir_i0i2 = 0. Since ffN^i = {r2 - ri)/V2 and pN-i = p 
we have aN~i — cti2 (compare eqs. (^) and (^9|)), so 
only derivatives with respect to ai2 remains. Thus it is 
convenient to introduce the notation 11^2 = ^'n-i- 

In the sum over angular wave function components in 
eq. ( |3^ ) only three different types of terms appear. As- 
suming 1 = 1 and J = 2 these types are classified by 
the set {fc, 1} either having two, one, or zero numbers co- 
inciding with the set {1,2}. Then eq. ( ^4|) is rewritten 
as 

0= (n22 + i;(ai2)- A)0(ai2)+ (35) 

N N N 



;=3 



/=3 



k>3,l>k 



with v{aij) — 2mp^V {y/2p sin aij) / . Multiplying this 
equation from the left by 0(ai2), followed by integration 
over all angular space except ai2 results in an integro- 
differential equation in a = ai2 of the form 

(n?2 + v{a) - a) 0(a) + v{a)2{N - 2) f dr ^(aia) 



+v{a)-{N ~2){N~3) J dr (^(034) = . (36) 

Here dr cx dr2jv_2 is the angular volume element exclud- 
ing the a-dependence; the normalization is J dr = 1. 
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This projection leaves for every value of a only two dif- 
ferent integrals due to symmetry between the first and 
second sum in eq. (^5|). Both of the remaining integrals 
can analytically be reduced to one dimension. The re- 
sults, collected in appendix H, are denoted by 



dT cf,{a3i) - Ri'^ 'V(a) , (37) 
J dT </.(ai3) = , (38) 

where '^^ is an operator acting on the function 4'{a) 

resulting in a new function of a. Mathematically R re- 
sembles a rotation operator, hence the choice of notation. 
Eq. (Bq) can now be written as 



0: 



+ - A + 2(iV - 2)v{a)RXi 



{N~2) 



+ -{N-2){N-i)v{a)R 



which is linear in the function 



{N-2) 
■34 



(39) 



3. Angular kinetic energy eigenfunctions 



We first consider eq. ( pQ]) for non-interacting particles, 
i.e. V = 0. Using the transformation in eq. ( p^ ) to get rid 
of first derivatives we find 



da 2 
9iV- 19 



V da^ 



(3iV-7)(37V-9) 



tan^ a 



- A U(a) = 



(40) 



where 
tion 



(a) is defined as the reduced angular wave func- 



(a) = sinacos(3^-^)/2^, 



(«), 



(41) 



in perfect analogy to the transformation from radial to 
reduced radial wave function for the two-body problem. 
Since <f>{a) for a physical state cannot diverge at a = or 
a = 7r/2, the boundary condition for the reduced angular 
wave function is ^(0) — 4){tt/2) = 0. 

The (non-reduced) solutions to eq. ( pO| ) is given by the 
Jacobi polynomials ||3^ 



_ p(l/2,(3W-8)/2) 



(cos 2a) , 



(42) 



where the hyperspherical quantum number K — 21/ = 
0, 2,4, . . . denotes the angular kinetic energy eigenfunc- 
tion with v nodes. The corresponding angular eigen- 
values are — K{K + 3A^ — 5). The lowest eigen- 
value is zero corresponding to a constant eigenfunction 



,(l/2,(3Af-8)/2) 




1. 



In fig. |i|a are shown examples of the reduced angular 
kinetic energy eigenfunctions for iV = 100 and the lowest 



three eigenvalues. The constant wave function, 4>k=0j is 
in the figure represented by (j)o (a) = sin a cos 



(37V-7)/2 



a. 



where |0o| then simply is the volume element in a-space. 
The oscillations are located at relatively small a- values. 
As seen in fig. ||b their amplitudes decrease as 
due to the centrifugal barrier proportional to tan^ a in 
eq. (|40|). Thus, as N increases the probability becomes 
increasingly concentrated in a smaller and smaller region 
of a-space around a = 0. 
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FIG. 1: The reduced angular wave function defined in 
eqs. @ and @, for a) = 100 and K = 0,2,4 and 
b) if = and = 10,100,1000. The normalization is 

/;/'da|0K(a)P = l. 

Some solutions may be spurious, i.e. each component 
is non-vanishing but the full wave function $ in eq. ( [3l| ) 
is identically zero. From eqs. (^) and (|3^) we see that 
the component for a spurious solution is an eigenfunc- 
tion of the angular kinetic energy operator. This special 
situation occurs for the K = 2 eigenfunction in eq. (|42" 
which satisfies 



N 

dr ^ <pK=2iaij) = . 

i<j 



(43) 



Solutions like 4'k=2 obtained by solving the full equation 
in eq. ( ^3| ) are therefore independent of the interactions 
and the eigenvalue is independent of p. 



E. The A-spectrum for large p 

For large values of p the short-range two-body poten- 
tial V with range h is non- vanishing only when a is smaller 
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than a few times h/ p. For larger values of a the "rota- 
tion" terms, Rcj), in the angular equation in eq. ( ^9| ) can 
then be omitted. 

For small values of a we find by substitution of r = 
-\/2psina instead of a in eq. (Pg) that the remaining 
equation (without rotation terms) can be rewritten as 



+ V{r)-E\u{r) , 



(44) 



where 2mp^E/h^ A + 9A^/2 - 9, u(^/2psina) = <^(a), 
and the reduced mass ^ — ■m/2. 

Each two-body bound state solution E <0 then corre- 
sponds to an eigenvalue A diverging towards — oo as —p^- 
Such solutions do not produce significant rotation terms, 
which is consistent with the omission in the derivation. 
The structure of the A^-body system is given by the fully 
symmetrized wave function of two particles in the bound 
state, while all other particles are far away thus produc- 
ing the large p. 

All other solutions to eq. (39) without bound two- 



body states correspond to wave functions distributed over 
larger regions of a-space. As the potentials then vanish 
for large p we are left with the free solutions, i.e. the 
free spectrum of non-negative A-values is obtained in this 
limit of large p in addition to the discussed diverging 
eigenvalues for bound two-body states. 

When the attractive potential contains precisely one 
bound two-body state with zero energy, the eigenvalue 
A approaches a negative constant for large p. This can 
be understood by considering a two-body state with en- 
ergy slightly below zero, which forces A to diverge slowly 
as — /O^. On the other hand, if the two-body system is 
slightly unbound, A instead converges to zero, which is 
the lowest eigenvalue of the free solutions. Precisely at 
the threshold it seems that A cannot decide and therefore 
remains constant. Thus for infinite two-body s-wave scat- 
tering length one A-value approaches a negative constant 
reached for large p, when the average distance between 
particles is much larger than the range of the interaction. 



III. ANGULAR VARIATIONAL EQUATION 

The angular equation is essential and solutions are not 
easily obtained. Proceeding with the Faddeev equations 
is one option, but we prefer first to derive the optimal 
angular equation within the Hilbert space defined by the 
form of the angular wave function in eq. (^l]). We also 
want to exploit that the two-body interaction is of very 
short range compared with the size of the system. 



A. Variational approach 

The angular Schrodinger equation for fixed p in eq. ( |2^ ) 
and the ansatz of the wave function in eq. (K^) allow us 



to express the eigenvalue as an expectation value, i.e 



A = 



1$) 



(45) 



Since is invariant with respect to interchange of par- 
ticles, the terms (0y |/in|$)j2 = {4>ki\hQ.\^) n are identical 
and eq. (BS) simplifies to 



»12|"i2 



'ki/n 



>ki n 



(46) 



The total angular volume element is dVlN-i = 
dr2i^~^^dri^^~^^cir27v-2j see eq. The integrands are 



independent of fl^l^ allowing to omit dfl^^ from 
the integrations. Using eq. (H) we then obtain 

N 

/ ^12 / dnN-2{hn - A) ^ (/.fci = . (47) 

'' k<l 

The wave function component (/)^2 is now varied until 
the lowest eigenvalue is obtained. This gives the integro- 
differential equation 



„ N r N 

/ d^lN-2 ^ (A^_i - \)4>ki + Vki ^ 

^ I. ^1 L. ™ 



k<l 



where the unknown functions, 



,(48) 



{otij), in fact 



all are the same identical functions of the different 



coordinates a; 



This result resembles the Faddeev- 



like equations of eq. (p3). Many terms are identical, 
e.g. / d^N-2 vi24>34 =~] dQN-2 vi24)56, since particles 
1 and 2 cannot di stin guish between other pairs of parti- 
cles, see appendix C 1 for the details. Collecting all terms 

yields 



J dnN-2 (n?2 + «i2 - a) (1)12 + G(r, aia) 



0,(49) 



where r denotes angular coordinates apart from ai2, and 
G{t, ai2) ^ ^712(1134 + v{ai2) + v{a3i) - A)(/)(q;34) 



'-n2v{a3i)<j){ai2) 



2niv{ai3){(f>{ai2) + (/'(a23)) 
+2r7,i(n^3 + v{ai2) + viais) - X)(f>{ai3) 
-|-n3(^-y(a34)((/)(a35) + ^(ais)) -I- w(ai3)0(a45)) 
+2n2v{ai3){(j){ai4) + 0(a24) + (f>{a34)) 
+2n2v{a34:)(l){ai3) + ^?^4w(a34)<?^'(a56) , (50) 

where ~ Hf^il^ ^ J ^ l)i and nf is defined from 



eq. (15) with k — N ~ I and with replaced by aij. 
In eq. (BG) all terms depend at most on coordinates of 



the six particles 1-6. The first three terms in eq. (49) do 
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not depend on the integration variables r leaving only 
G'(r, ai2) for integration. 

By appropriate choices |0 of Jacobi systems, the 
relevant degrees of freedom can be expressed in terms 
of the five vectors fjj^^i, . . . ,fj]S!-z. One is the argu- 
ment of the variational function and not an integration 
variable. The remaining twelve-dimensional integral is 
then evaluated with the corresponding volume element 
dr oc n?=2 '^^"^ dilif^ where the normalization is 
J dr ~1. Then eq. ( ^ ) becomes 



(51) 



where we used that the first terms are independent of 
the integration variables. Eq. ( |5l| ) is a linear integro- 
differential equation in one variable cont aining up to fivc- 
dimensional integrals, see appendix |C 2|. 



B. Simple models for short-range interactions 

The short-range two-body interaction with s-wave 
scattering length Os has in mean-field contexts Q been 
modelled by the three-dimensional (5-function potential 



Vim) 



m 



(52) 



With the constant angular wave function, <f>x=o ~ 
'^i'^j <j>K=oioiij), the expectation value of the angular 
Hamiltonian /iji becomes 



Xk=o = {^K=o\hn\^K=Q)n 

N 
k<l 



(53) 



without contribution from angular kinetic energy. With 
the (5- interaction in eq. ( ^2|) we obtain (see also |20t] ) 



A 



K=0 



2 r 



TT r 



N{N -1) 



(54) 



The (5-interaction however does not contain the pos- 
sibility of studying the short-range properties such as 
bound two-body systems and collapse of a condensate 
into clusters. An improvement is made by using a finite- 
range, but still short-range, potential. When p is much 
larger than the potential range, &, we find an expectation 
value of the angular energy of the same form as A^^q in 
eq. (HI): 



A 



finite P»6 
K=0 * 



-r(^) 



(55) 



The strength is now collected in the parameter ob instead 
of a=, where 



Ob 



d^Tki V{rki) 



(56) 



is the Born-approximation to the scattering length as- 

However, in the opposite limit, when p <^ b, the re- 
sult is strongly dependent on the shape of the poten- 
tial. For example, for a Gaussian potential V(rki) = 
Vb exp(— r^;/5^) with as = ^/TTmb^Vo/iAfi^) we obtain 



\ finite P<'' 



As seen from these two limits there are some interesting 
scaling-properties for finite-range potentials. The angu- 
lar eigenvalue at a given iV- value depends only on a^/b 
and p/b. More specifically we find for a Gaussian poten- 
tial 



Vkl 



y^b\b 



3-2(p/&) sin" Qfci 



(58) 



which implies that for a given a^/b the angular eigenvalue 
A is a function oi p/b only. Moreover the radial potential 
can be written as 



2rab^U 



A (3A^-4)(37V-6) (p/6)2 



A{p/bf 



(59) 



and the scaled energy, 2mb'^E/h^, is then for a given iV- 
value a function of ob/^ and bt/b only. These scaling 
properties are useful in model calculations. 

The scattering lengths are apparently essential param- 
eters. We show in fig. ^ the s-wave scattering length, Og, 
as a function of the strength parameter |aB|/& for both 
the Gaussian potential and for an attractive box poten- 
tial. When the numerical value of the parameter as is 
increased from zero, the scattering length varies slowly 
and roughly linearly with ob for small ob until the value 
where diverges as a signal of the appearance of 
the first two-body bound state. The threshold value of 
ag ^ is different for the Gaussian and for the box poten- 
tials, but flg/aB as a function of OB/og'' results in virtu- 
ally the same curves, see fig. |^b. This indicates that for 
simple potentials, the behaviour is approximately shape 
independent. 



C. Short-range approximation 

The angular eigenvalue equation, eq. (|T]) , simplifies in 
the limit when the two-body interaction range b is much 
smaller than p. Then the integrals are either analytical 
or reduce to one-dimensional integrals. This substitution 
of a (5-function is only allowed for the potentials appear- 
ing under the integrals. Otherwise the Thomas collapse 
occurs Thus apart from the local terms containing 
v{a), the results only depend on ob defined in eq. ([s^). 
For example the J dr w(Q;34)-term reduces to 



dr v{a34) ~ vi{a) 



2 r( 



r(^ 



cos'^ a 



(60) 
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FIG. 2: a) Scattering length as divided by the typical po- 
tential range 6 as a lunction ol divided by b. b) Scatter- 
ing length as divided by as a function of divided by 
where the first bound state occurs. Results are shown 



for the Gaussian potential V(r) = Vbexp(— r /& ), (ajj /fe = 
-1.1893), and for the box potential V{r) = VoQ(r < b), 
(a^'Vfc = -0.82248). 



Similarly the / dr u(ai3)-terni reduces to 



dr v{ai3) — V2ia) = 



3^3 



(61) 



where sin/3o = tana/v^ and 8 is the truth function. 
The remaining terms can in this Umit be expressed 
through vi{a), V2{a), R^^'^ from eqs. (| 



other related rotation operators R) 



definitions are given in appendix C 3 



(n) 
kl- 



) and and 
Corresponding 



The reductions can be understood qualitatively via 
the diagrams in fig. |3[ which shows the geometry when 
the zero-range interaction is contributing to the inte- 
grals. As an example, in the integral j dr v{aiz)4>(a3A), 
see fig. particles 1 and 3 must be close together 
as shown fig. ^d. Then the distance between parti- 
cles 3 and 4 appearing in ^34 is approximately equal 
to the distance between particles 1 and 4. Therefore 
/ dr w(ai3)(/'(a34) ^ J dr v{ai3)(t){ai4). 





P34 ~ <Pl4 
Vl3 ~ S{f 13) 

a) 3b)'" 
FIG. 3: Simplifications due to short-range potentials. 

For & <C p we get almost independent of the potential 
dr G{T,a) ~ (^-YVi{a) + AniV2{a)j (l){a) 



34 



v(l){a) + 2niR[^ ^^v(j){a) 



+ Y(i?^r'^ni40(a) + {v{a) A)i?(r^V(a) 



-|-2riirR 



£,(A'-2)jj2 



^•13 



i3V\a) + iv{a) - A)i?i3 



{N-2) 



0(a) 



(JV-3)^ 



n-4 

-f— wi(a)i?34 
+n3V2{a)R[^;^^^(f){a) 



(a) + n3Vi(a)-R; 



3435 'i 



(62) 



+n3Vi{a)R[^ </>(") + 2n2Ui(a)^^4\3(/)(a) 
+2n2V2{a)[2R^^^,^<j,{a) + Rg\^(b{a)] . 

In the extreme zero-range limit the two terms R^ f '^'^v(f> 
are proportional to 0(0). Using eqs. ( |60| ) and (61) we 
define v and g by 

u(a) EE ^i;i(a) , g{a)v{a) = 2niV2{a) (63) 



and rewrite eq. (|62|) as 

dr G{t, a) ~ v{a) (1 + 25(a)) 0(a) 4- ^0(a) 

+i)(a)(l + 5(«))0(O) , (64) 

where R(j){a) collectively denotes the remaining rota- 
tional terms. Collecting the terms from eq. ( |5l| ) then 
leads to an integro-differential equation in one variable: 



= (n?2 + via)+ «(a) (1 -I- 25(a)) - a) 0(a) 

+7?0(a)-f {)(a)(l + 5(a)) 0(0) . (65) 



Eq. (|6g) is a linear eigenvalue equation (in 0) for one 
variable a. The eigenvalue A(p) is the key quantity in 
the much simpler radial equation. 

If we neglect the rotational terms in R, assume 0(0) ~ 
0, and use that 5(a) <C 1 when N is large, we get a so- 
lution A « v{0). Since v{0) = Affi*^ at large N, this 
eigenvalue is the same as found in eq. (|5^) from the ex- 
pectation value of a finite range potential in a constant 
angular wave function. Eq. (65) thus predicts an eigen- 
value proportional to p^^ with a proportionality factor 
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given by eq. (p5|), where the determining parameter is ob- 
However, as discussed in [|o) the angular eigenvalue for 
small scattering length as is determined by as only. This 
is for a sufficiently weak interaction the correct limit for 
Ob as discussed in connection with fig. ^. 

The repulsive effective potential in w(a) pushes the 
wave function into a narrow region at small a outside 
the repulsive core of w(a). Such a solution is often not 
a good approximation because the rotational terms are 
important. Also for attractive potentials the confinement 
to small a can not be achieved by the v{a) term. It is 
therefore crucial to include all the terms of eq. (|65| ) as 
discussed in [^. 



IV. NUMERICAL ILLUSTRATIONS 

The method has to be tested by solving the derived 
equations for realistic parameter choices. We first discuss 
qualitatively which numerical technique we apply. As we 
shall explain the difficulties increase as N and p increase. 
We shall therefore concentrate on relatively small values 
of N, which happens to be a region of growing inter- 
est in the art of making Bose-Einstein condensates. The 
first physical results focus on the angular eigenvalues and 
the related wave functions. Then the fully defined radial 
equation and its solutions are discussed. 



A. Method of solving 

A usual method within the hyperspherical formalism 
is to expand the angular wave function on kinetic en- 
ergy eigenfunctions, the so-called hyperspherical harmon- 
ics [po|]3^. However, since the hyperspherical harmonics 
contain oscillations at about a ^ l/K, large iC's of the 
order of Kj^ax ^ p/b are needed to describe potentials 
limited to a < b/p. This i^max becomes very large for 
the application to Bose-Einstein condensation, so a basis 
of hyperspherical harmonics is not suitable in this con- 
text. Thus using only K = for a zero-range interaction 
|pO| must be far from the optimal solution. 

Instead we choose a basis of discrete mesh 
points distributed in a-space — > = 
[(/)(q!i ),..., 0(a„i)i to take into account 

the short range of the potential and to keep suffi- 
cient information about small a. Derivatives are then 
written as finite differences |39| and integrations like 



R<j){oi) of eq. (65) can be expressed in matrix form, 
i.e. R(t)(a) ^ Rcj). 

Numerical computation of the integrals becomes in- 
creasingly difficult with decreasing interaction range. 
This can be understood in terms of the ai2 coordinate, 
since the potential at a given p and a given range, b, 
of the interaction, is confined to an ai2-region of size 
Aai2 ~ b/p, which for Bose-Einstein condensates easily 
becomes very small and thus cannot be handled directly 
numerically. In the following we shall use the equations 



obtained in the short-range approximation, where the dif- 
ficulties are much smaller and the physics content is still 
maintained. 



B. Angular solutions 

The behaviour of the lowest eigenvalues A depend 
strongly on the potential. The characteristic feature is 
the large distance asymptotic behaviour, i.e. divergence 

as 



eq. 



^ corresponding to a bound two-body state (see 
4)) or convergence as p~^ towards a finite value 



corresponding to the spectrum for free particles. The 
constant of proportionality to p~^ is qualitatively recov- 
ered as the predicted dependence on Og. The excep- 
tion arises for infinite scattering length at the threshold 
for two-body binding, where one angular eigenvalue at 
large p approaches a negative constant. This structure is 
illustrated in fig. ^, where the lowest angular eigenvalue is 
shown for various Gaussian potential strengths covering 
the region from unbound to bound two-body states. 
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FIG. 4: Angular eigenvalues for A^ = 20 and various parame- 
ters {a-g/b, as/b) as shown on the figure. A star refers to the 
first excited state. 

For repulsive potentials the eigenvalues are positive 
and the lowest approaches zero from above. For weak at- 
tractions the lowest A is negative approaching zero from 
below. When the attraction can bind a pair of particles 
the divergence is seen while the second A for the same 
potential is positive and approaching zero from above in 
qualitative agreement with the lowest eigenvalue for a 
repulsive potential. The higher eigenvalues would then 
converge to K{K + 3N ~ 5) as 1/p, where if = 4, 6, 8.... 
The solution ioi K — 2 is not allowed, corresponding to 
removal of the non-physical spurious solution of the Fad- 
deev equations. Increasing the attraction to allow two 
bound two-body states would then shift the asymptotic 
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spectrum such that one more eigenvalue diverge while 
the non-negative energy spectrum remains unchanged. 

At each threshold for the appearance of a new bound 
two-body state one eigenvalue asymptotically approaches 
a negative constant. This especially interesting eigen- 
value is responsible for the structure of the A^-body sys- 
tem for very large scattering lengths. 

The total angular wave functions are determined as a 
sum of the components. We show in fig. || an example 
of the reduced wave function for a potential with one 
bound two-body state. The amplitude increases with 
p and concentrates at very small values of a. This re- 
flects the convergence towards the two-body bound state 
in agreement with the transformation ri2 = \/2psina. 
Recovering this behaviour numerically is essential, oth- 
erwise the large distance behaviour cannot be described. 
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FIG. 5: The lowest angular wave functions for A = 20 and 
aB = —1.506, Os = 5.316 for three values of the hyperradius. 
This potential has one bound two-body state. 

The angular eigenfunction varies with the strength of 
the interaction. Examples of this variation are shown in 
fig. 1^. The non-interacting wave function has only nodes 
at the cndpoints. The repulsive case shows an oscillation, 
which lowers the angular energy due to the rotation terms 
PH . The fast change at small a is typical for interacting 
particles. The wave function for the excited state has 
an additional node. The corresponding lower-lying wave 
function is shown in fig. ^ 

The wave function for infinite scattering length corre- 
sponds to an interaction where the two-body bound state 
is at the threshold for occurrence. This eigenfunction re- 
sembles those where a bound two-body state is present, 
compare with the results shown in fig. |^ with different 
scales on both axes. However, now (thick curve of fig. |^) 
the wave function is located at larger a- values and a node 
is present in the tail at an intermediate a ^ 0.25. 

The structure of the component of the angular wave 
function is further illustrated by the second moment de- 
fined by (^12)0 = 2p^(</)| sin^ al^). A number of these 
moments for different interactions are shown in fig. ^ as 
functions of p. For states obtained from repulsive po- 
tentials, moderately attractive potentials without bound 
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FIG. 6: Angular wave functions for A = 20 and p = 5006 
for different interaction parameters (03/6,03/6) as shown on 
the figure. The A = curve corresponds to a non-interacting 
system. A star refers to the first excited state. 

two-body states, and for excited states of positive A, this 
moment, (r^j)^: increases proportional to p^ for large p. 
This resembles the behaviour of the expectation value 
in the lowest angular state for a non-interacting system, 
a: = 0, where (rJa)^ = '^P^/i^ - !)■ The qualitative 
explanation is that large p implies the limit of a non- 
interacting spectrum with the corresponding wave func- 
tions. The particles exploit the possibility of being far 
from each other, since there is no energetic advantage of 
being close due to the lack of a bound two-body state. 

In contrast a different behaviour is observed when the 
potential can bind two particles, i.e. (r^j)^ approaches 
a constant at large p. This can be understood as the 
angular equation in this limit approaches the two-body 
equation in eq. (^). The radial wave function in the 
zero-range limit converges to u{r) = e~^^°'^, where as is 
the scattering length. The second moment is then found 
as (ulr^ju) — a1/2, which in the limit of large p reproduce 
(''12)0 for ^s/b — 9.32 and ag/h = 5.31, see fig. 0. 

This can be expressed in a different way: when a two- 
body bound state is present, the angular wave function 
is at increasing p squeezed inside the potential, since the 
range in a-space decreases proportional to p^^, and we 
obtain (0| sin^ Ci\4') (x 1/p^. The distance between a pair 
of particles is therefore independent of p at large values 
of p. This means that pairwise the two-body bound state 
is approached while all other particles are far away. The 
symmetrization does not affect this conclusion. 

At the threshold for two-body binding, infinite scatter- 
ing length, we again observe the intermediate behaviour 
resembling a logarithmic dependence in fig. |^. 



C. Radial solutions 

The most interesting angular eigenvalues either con- 
verge to zero or remain constant for large p. We therefore 
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FIG. 7: The second moment (r\2)(t> as a function of hyper- 
radius for A'^ = 20 for solutions to the angular variational 
equation with different interaction parameters specified in the 
figure by (ag/fe, as/b). Also shown is the K = 0- value. A star 
refers to the first excited state. 

select an interaction where A approaches zero relatively 
slowly from below. This is then a weakly attractive po- 
tential without bound two-body states although not very 
far from this threshold. The parameter as/h = —15.7 cor- 
responds to this case. The resulting A shown in fig. ^ is 
used to compute the radial potential in eq. ([2^). The 
total potential is shown in fig. ||. 

The radial potential always diverges at large p due to 
the harmonic external field and at small p due to the 
centrifugal barrier term. Thus there is always infinitely 
many bound states. For the moderate attraction used 
in fig. H the potential for the lowest angular eigenvalue 
has a global negative minimum at small p separated by a 
barrier at intermediate p from a local positive minimum 
at p - pt = bf The minimum at small p disap- 

pears quickly with a slightly less attractive potential. It 
becomes deeper and wider for a more attractive poten- 
tial, where the barrier at intermediate p simultaneously 
is reduced and eventually disappears completely. 

The radial potential corresponding to the second adia- 
batic angular potential is also shown in fig. ||. Since this 
A is positive for all p an attractive pocket {U < 0) in the 
radial potential cannot arise. It coincides with the lowest 
radial potential for large p and due to the lack of attrac- 
tion at small p the potential therefore diverges to +oo 
for small p without going through another minimum. 

The coupling terms of eq. ( p5| ) contribute at most 
about 1 % compared with other terms of the full radial 
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FIG. 8: Radial potentials Uo and Ui from eq. (g6|) corre- 
sponding to the two lowest angular potentials for = 20 
and as/6 = —15.7 in fig. ^. We model the experimentally 
studied systems [j40| of *^Rb-atoms with oscillator frequency 
1/ = uj/{2-k) = 205 Hz and interaction range 6 = 10 a.u., thus 
yielding bt = \Jhl(m,uj) = 14426. Also shown are the two 
lowest energies and the radial eigenfunctions /gs and /cx for 
the lowest radial potential in the uncoupled radial equation, 
eq. d^). 

equation, eq. (p^, and can be omitted. We therefore ar- 
rive at solving the uncoupled radial equation, eq. (|2^), 
the solutions of which are considered in the following. 

The external field is negligible when p <^ VNbt and 
the radial potential is therefore negative when A + {3N — 
4)(3iV — 6)/4 < and p is sufficiently small. Then gen- 
uinely bound many-body states of negative energy are 
possible in our model without the influence of confine- 
ment from the trap. The radial equation corresponding 
to the parameters of the lowest potential shown in fig. g 
has only one negative-energy solution with the wave func- 
tion located in the global negative minimum. 

The first of the infinitely many excited states in this 
potential is located in the local minimum at larger p cre- 
ated by competition between the centrifugal barrier and 
the external harmonic oscillator potential. This excited 
state is usually referred to as the Bose-Einstein conden- 
sate, because the corresponding wave function is similar 
to that obtained in experiments (2| . It also resembles the 
wave function found from the purely repulsive radial po- 
tentials arising from both the second A with the present 
attractive interaction and from repulsive two-body inter- 
actions. 

From eq. (|) we have (p^) = N{rf)~N{R'^) = N{rf) - 
ihl/2. Using eq. (|) we get 2(p2) = (TV - l){rl^). The 
root mean square radius {rf)^^"^ of the condensate state, 
/ox, is then 2 % larger and the energy E'cx is 3 % lower 
than the Gross-Pitaevskii results. The lowest radial state 
in Ui has both energy and root mean square radius about 
10 % larger than the Gross-Pitaevskii results. 

The moderate attraction used to obtain the potential 
in fig. ^ is therefore seen to produce a lower-lying many- 
body bound state with an average distance between the 
particles about 37 times smaller than in the condensate 
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state. The structure of this state could as well be char- 
acterized as a condensate (condensed A^-body state), but 
it is much more unstable due to the many orders of mag- 
nitude larger density and the subsequent larger recom- 
bination probability This ground state, /gs, has no 
parallel in usual mean-field computations. 

Even the state with Egs < 0, (p^)^/^ ~ 2166, has a root 
mean square distance between two particles, (''12)^^^ ~ 
70&, much larger than the interaction range b. This is 
a sufficiently large average distance between particles to 
assume that the short range details of the two-body po- 
tentials are unimportant. Increasing the attraction more 
negative-energy states appear in the attractive pocket in- 
side the external trap. They occur at increasingly larger 
densities, but still low enough for the short-range details 
to be unimportant, see more details in [ p^ . 

At the two-body threshold, when the s-wave scatter- 
ing length diverges, potentially infinitely many negative- 
energy many-body states could occur. However, the ex- 
ternal harmonic oscillator field limits the number iV^; of 
these new negative-energy states, located inside the trap 
and outside the two-body potential, to approximately 

Ne » O.STVln (J^^ , for h « iV|as| . (66) 

The number of these multi-particle Efimov states there- 
fore scales linearly with N and logarithmically with the 
ratio of trap length and interaction range, see also |T2[ . 

The stability of the A^-body system as a coherent ob- 
ject (condensate) depends on the radial potential. Start- 
ing with a condensate in the minimum of /ox in fig. ||, 
decay occurs both through two- and three-body recom- 
bination processes and through macroscopical tunnel- 
ing through the barrier to state(s) in the global mini- 
mum. These denser states recombine faster due to the 
smaller interparticle distances. Increasing the attraction, 
the rate for macroscopical tunneling increases due to a 
smaller barrier. Eventually, at large (infinite) scattering 
length, the barrier has vanished and the A^-body system 
can contract freely, i.e. populate the many-body Efimov 
states. The time scale for this contraction is estimated 
to be about 0.1 ms |p^ , which can be compared with ex- 
periments [ pT| and mean-field computations pO[ |. Thus 
sudden removal of the barrier over the first short period 
of time lead to spreading of the probability to the smaller 
distances where recombination then takes place and the 
condensate decays with a corresponding lifetime. This is 
qualitatively in complete agreement with the measured 
decay function (ll|. 

V. SUMMARY AND CONCLUSIONS 

We have formulated in details a new method to inves- 
tigate two-body correlations within symmetric A''-body 
systems. In contrast to the product assumption for the 
wave function in the mean-field approximation we use a 
Faddeev-type decomposition. The overcomplete basis is 



simplified by using in each of the two-body components 
only the lowest possible number of partial waves, i.e. s- 
waves. The allowed Hilbert space is then dramatically 
reduced, it is not complete, but hopefully precisely de- 
signed to treat the two-body correlations. 

We use an adiabatic hyperspherical expansion with 
the hyperradius as the adiabatic coordinate. We de- 
rive the optimal angular equation for the two-body am- 
plitude (the Faddeev component) arriving at an one- 
dimensional angular integro-differential equation. Its 
eigenvalues are closely related to the effective hyperra- 
dial potential, which receives contributions from the two- 
body interactions, the kinetic energy operator and from 
the external field confining the system. This potential 
diverges both at small and large hyperradii due to the 
centrifugal barrier and the external trap, respectively. 

For repulsive two-body interactions only one minimum 
in the radial potential exists, but for moderately attrac- 
tive two-body potentials two minima are present sepa- 
rated by a barrier. The minimum at shortest hyperradii 
is at negative potential and is therefore able to bind the 
system without help from the confining external field. 
For a state in this minimum the particles are on aver- 
age still far outside the range of interaction with each 
other. The other minimum at larger hyperradii can also 
be pronounced enough to hold localized stationary states, 
which are similar to the solutions usually referred to as 
Bose-Einstein condensates. 

As the attraction increases and approaches the thresh- 
old for two-body binding, i.e. the scattering length in- 
creases towards infinity, the intermediate barrier disap- 
pears and only the negative minimum at smaller hyper- 
radii remains. However, it becomes both deeper and 
wider and as a consequence more bound states of negative 
energy appear. These many-body states have character- 
istic features of Efimov states. They have inter-particle 
distances much larger than the range of the two-body in- 
teraction and are therefore universal structures indepen- 
dent of details of the potential. Their number increases 
proportional to the number of particles and logarithmi- 
cally with the ratio of trap length to interaction range. 

In conclusion, we have formulated a method to treat 
two-body correlations in an iV-body system. We ap- 
plied the method to Bose-Einstein condensates and ob- 
tained a simple one-dimensional pictorial description in 
terms of an effective length coordinate. Macroscopic col- 
lapse is conjectured to proceed via the new universal Efi- 
mov states at intermediate hyperradii, which quickly, due 
to the much larger density, subsequently recombine into 
dimer or trimer states. This decay can be studied quan- 
titatively in the present model. Another unique feature 
of the model is to provide a solution even in the case of 
very large scattering lengths, where the Gross-Pitaevskii 
equation breaks. 
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APPENDIX A: COORDINATES 



fined by 



For use in the calculation of matrix elements different 
Jacobi trees can be chosen . The relevant ones in this 
context are shown in fig. ^ 



VN-i = ■^(f^2 -ri) , 



■nN-2 = 



/f(r3-i(f2+rl)) 



m-3 = ^ ^{fi - ^{fs + f2 + fl) ) , 



(Ala) 
(Alb) 
(Ale) 



Vn-i 



a) 




VN-3 



FIG. 9: Jacobi trees: a) standard, b) (12)(34), c) (123)(45), 
d) (12)(345), and e) (12)(34)(56). 



The coordinates of the standard tree of fig. m are de- 



m = 



N -1 



{rN-i + ■■■ + ri)) - (Aid) 



In the (12)(34)-tree of fig. two of the vectors are 
different from the standard tree: 



rfN-2 = -^(ri-ra) , (A2a) 
?7Ar_3 = i(f4 + rg - - ^i) • (A2b) 

In the (123)(45)-tree of fig. ^ two of the vectors differ 
from the standard tree: 



m-3 = -^(rs - ^4) 



VN-i 



(A3a) 

(^(^5 + - ^{rs +r2 + n)) . (A3b) 



6/1 



In the (12)(345)-tree of fig. ^ three of the vectors de- 
viate from the standard tree: 



m-2 = ;^(^4 - ?%) , 



(A4a) 
(A4b) 



-4 = ■\/f (^(^5 +ri + fa) - i(f2 + fl)) . (A4c) 



In the (12)(34)(56)-tree of fig. ^ four of the vectors 
are different from the standard tree: 



m-2 = -^(f4 - fa) , TIN -3 

VN-i ^ -^(ri + r3 - r2 ~ r2) , 



V2 



(fg-fs) , (A5a) 



4 / f 6 + f 5 f 4 + f 3 + f 2 + f 1 



(A5b) 

(A5c) 



Since only inter-relations between ttat-i, ?fAr-2, and 
f]N-3 are needed in evaluating the matrix elements, we 
will use the common notation: 

ryAr_i = psina , 7yjv-2 = pcosasin/3 (A6) 
W-3 = pcosacos/3sin7 , fjk ■ m ^ rjum cos dk,i , (A7) 

where Ok^i is the angle between the fc'th and Tth Jacobi 
vectors. We abbreviate 6'Ar_i^Ar_2 ^x, Qn-i.n-3 ~> (^y, 
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and 0N-2,N-3 dz- The azimuthal angle Lp determining 
the projection of r/jv-s onto the plane of ^^v-i and r/Ar_2 
is defined such that 

cos 9z = sin 9x sin 6y cos </? + cos 9x cos 9y . (A8) 



With r = 7, 6*2:, 6*^, (/?} a matrix element of an arbi- 
trary function / of all the variables a and r then becomes 



dr f{a,T) 



fj' dl3 sin2 f3 cos3^-io p fj' d-f sin^ 7 cos^^-^^ 7 sin 0, /J^ d^,, sin 0^ d^ f{a, 

i9yjl'" dip 



/o"/' d/3 sin^ /3 cos3^-io /3 /g"/" d-f sin^ 7 cosS^"" 7 d9x sin 0, /p" d9y 



n7r/2 



,(A9) 



where the normalization is J dr = 1. We need relations 
for interparticle distances and define ffij = (r 
and the angle a^- related to r]ij = p sin aij = 1 



APPENDIX B: INTEGRALS (FADDEEV) 



where the last notation is convenient for repeated use. 

To describe three particles in /c?t ^(aia) simultane- 
ously, Jacobi vectors of the standard tree are needed. 
The distance between particles 1 and 3 is related to the 
corresponding Jacobi vector 



Eqs. ( pTj) and (|3|) are evaluated as follows. 

In the integral / dr 0(034) a convenient choice of co- 
ordinates is th e al ternative Jacobi (12)(34)-tree of fig. ^ 
given by eqs. ( |A2[ ). The angle 034 is associated with the 
distance r34 — \/2rj34 by the relation 

?734 = VN-2 — p cos a sin /3 = p sin 034 (Bl) 
sin 034 = cos a sin /3 . (B2) 

The integrand, 0(034), only depends on 034 , w hich is a 
function of a and f3. Therefore at fixed a eq. ( A9) reduces 
to 



dr 0(034) 



4 r(^) 



J^^' dp sin^ p cos3^-io P 0(034) 
Jo^^dp sin2/3cos3^'-iO/3 



/^r(^)7o 



7r/2 



dp sin /3cos' 



3Af-10 



P <?!>(a34) 



R 



{N~2) 
'34 



0(0) 



(B3) 



Vl3 



V2 



1 \/3 

(r's - n) ttW-i + , (B4) 



The hyperangle 013, associated with the distance be- 
tween particles 1 and 3, through 7713 — r 13/^/2 — 
psinoi3, is then 



sm oi3 



■ sin'^ o 



■ cos o sin P + 



(B5) 



sin a cos o sin P cos 9^ , 
where 9x is the angle between the Jacobi vectors 77Ar_i 
and ffpf-2- Note that 0(013), through 013, for fixed o 
depends on P and 9^, which leaves a two-dimensional 
integral. Therefore eq. ([Ag|) becomes 



/ 



J^^'^dp sin^ Pcos^^-^° PJ^ d9x sin 61^ 0(013) 



dr 0(013) ^ 



J^^^ dp sin^ p cos3^-io /? /J" d6l:r sin 9^ 



2 r(^) 



dp sin^/3cos3^-i°/3 / ^6*^ sin 



iai3j 



(B6) 



This integral can be reduced to one dimension by a partial integration. The final one-dimensional integral becomes 



dr 0(oi 



4 r 



2 ^ sin-iocos«-3^o 



•3Af~7 



7r/2— |7r/6— a 



(ioi3 COS 7 sinoi3cosai30(ai3) 



(Q-7r/3)e(Q>7r/3) 



(7r/3-a)e(7r/3>a) 



doi3 cos'^^ ^7 sinoi3 cosai30(ai3) 



(B7) 
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where sin^ 7^ — 4(sin^ a + sin^ ai3=FsinasinQ!i3)/3, and 
Q is the truth function. The last notation is convenient 
for repeated use. The integration regions of the integrals 
are shown in fig. ^ At fixed a the range of ai3 in the 
first integral is over regions I and II, while in the second 
it is only over region II. For = 3 the integrands of the 
two integrals are identical and thus the integration over 
region II in the first integral cancels the second integral. 




FIG. 10: The regions of integration in eq. (B7) 



APPENDIX C: INTEGRALS (VARIATIONAL) 



We first divide the integrals of eq. (48) into similar 
terms, then compute them in general, and finally in the 
zero-range limit. 



1. Counting terms 

We have to evaluate the double sums of eq. (^) in- 
cluding the potential: 



N 



N 



(Cl) 



k<l 



Three types of terms occur, due to the fact that we vary 
the wave function component 0*2 in eq. (^): the poten- 
tial concerning particles 1 and 2, the potential concerning 
one of the particles 1 or 2 and a third particle and the 
potential concerning neither particle 1 nor 2, but a third 
and a fourth particle. We obtain 



N N N N 

'^Vkl ^ V12 +^Vu +^V2l + ^ Vkl 
k<l 1=3 1=3 k>3,l>k 

^ vi2 + 2(iV - 2)i;i3 + ^{N -2)iN~ 3)v3i , (C2) 

where the arrow indicates the identity of the terms after 
integration over all angles except ai2 (analogously to the 
steps leading up to eq. (|3^)). Treating each of these in 
the quadruple sum, where the repeated use of arrows (^) 
has the meaning given just above: 



Fixing 0*2 and V12 yields three different terms: 

N 

V12 ^ 4>nin = 
m<n 

N 

■"12(012 + 2^(Pln-^- / , 92 

n— 3 n— 3 m>3,n>m 



(C3) 



! (012 + X! "^l" + X! "^^n + ^ 
n=3 n=3 m>3,n: 

vi2 (012 + 2(iV ~ 2)013 + i(iV - 2)(7V - 3)034) 
as shown in fig. pT| . 



<j>*,V,(j> 



b) 



3 c) 1 



FIG. 11: Illustration of 12^12-terms. 

Fixing 0J2 and W13 yields seven different terms. These 
can be identified in two steps, the first of which separates 
into four different sums: 



JV 

vi3 < 

N 



(C4) 



Vl3 



N N N N 

( ^ 01n + ^ 02n + ^ 03ri + ^ 



n— 2 n— 3 n— 4 m>4,n>m 

Each of these four terms are then identified as: 



N 



Vl3 ^1 



n=2 



Wl3 U'12 



N 

h3 + ^ 01n 

n=4 



^13(012 + 013 + {N - 3)014) 



(C5) 



The similarity of the terms in the first sum becomes ap- 
parent when carrying out the integration over dQN~2, 
e.g. / dil7v-2i^i30i4 = / dflN-2Vi34>i5- The rest are 
found in similar ways: 



N N 
Vl3 X! '?^2n = Vi3 (023 + ^ ' 



n— 3 



Wl3 (023 + (A^- 3)024) , 



N 



Vl3 ^ 03n Vi3{N - 3)034 
n=4 

N 



(C6) 
(C7) 



E 



1 



Vl3 / ^ H'mrL 
m>4:.n>rn 



"i3 2(^-3)(Af- 4)045 .(C8) 



The resulting seven types are shown in fig. 

Fixing 012 and V34 yields six different terms, identified 
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See the six types in fig. O. 



2. Evaluation of terms 



g) ^ ^ 3 
FIG. 12: Illustration of (/)j2i'i3-terms. 

as follows. The first step is: 

N N 

■^34 ^ 0mn =■^^34^^ 4>ln + (C9) 
m<n n— 2 

N N N N 

02n + ^ '?!'3ri + ^ </>4n + ^ 0mn , 
n=3 71=4 n=5 m>5,n>m 

In the next step the sums are treated: 



The term of fig. |ll| a is trivial since the integrand is 
independent of r. The terms of figs. |ll|b |ll|c, p^a, |lT 



,4 2, 



• 

3 y 



3 b) 1 
.4 2, 



AT 

1'34 ^ (/iln = W34 + </'13 + '/'14 

n— 2 n— 5 



-> «34(012 + 2<?ii3 + (iV - 4)015) 
Af AT 

V34 ^2 = "^34 (^23 + '?!'24 + ^ 4>2n 

n— 3 Ti— 5 

->'y34(2?!)i3 + (iV-4)(/)i5) , 

^34 X! = ^^34 ^034 + ^ (y^3n 

n— 4 n— 5 

->'y34(034 + (Af- 4)035) , 
Af 

^34 X! 1^34 - 4)03 



^35 



(CIO) 

(Cll) 

(C12) 
(C13) 



3 d) 1 

4 2, 



e) ' 



3 1 3 
FIG. 13: Illustration of 0i2t^34-terms. 



|13|a, and |13| c can be evaluated by eqs. (|B3D and ( |B7| ) 



n— 5 

A? 

W34 ^ 0mn 

m>5,n>m 



1 



The term of fig. 12b becomes with the use of the stan- 



^^34 2 - 4)(iV - 5)056 • (C14) dard Jacobi tree of fig. |9|a 



dT /(ai3) 5(a23) = -r ^ 3jv% / sin' cos^^-^" /3 / sin /{au) 3(023) , 

v"" 1 ( — 2 — ) "'0 Jo 



(C15) 
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3 1 \/3 

sin^ ai(2) 3 ~ ^ ''^^^ ^^^^ ^ ^ ^^^^ ^ ~^ ~2 ^^'^ ^ ^'^^ ^ ' (^-'^^) 

The term of fig. p^ f becomes with the use of the alternative (12)(34)(56)-tree of fig. ^ 

dT w(a34) 0(a56) = / sin^ /Jeos^^-^^ /3 / ^7 sin^ 7 cos^^"" 7 ^(034) ^(ase) , (C17) 



TT 

sin a34 = cos a sin /3 , sin asg = cos a cos /3 sin 7 , Aat = (3iV — 8)(3A^ — 10)(37V — 12) . (C18) 
The terms of figs, and are evaluated using the (123)(45)- and (12)(345)-trees of figs, ^c and so 

-^N I ,^ . O ^ 'i AT 1 n / 



J dTh{a,T) = ^ dp sin^ Pcos^'^-^° P d-y sin^ 7 cos^^^^^ ^ ^ dO^^^ s\u9^^, h{a,' 



(C19) 



where I^{a,T) can be either v[azA)4'{az^) or /(ai3)5(a45). The relevant angles are sin 0:34 = cosasin/?, sin^ a35 = 
cos^ a(3 cos^ /3 sin^ 7 + sin^ /3 + 2-\/3cos/3sin/3sin7Cos6'2)/4, sina45 = cos a cos /3 sin 7, and ais given by eq. (p^). 
Note that the integral / dr v{a^4)4>(pii^) = J dr (/)(ai3) v{a45). 

The terms of figs, p^, Oe, p^, and Oe are evaluated using the standard Jacobi tree. Then eq. ([A9|) reduces to 



dr f{ai3) 5(az4) = 74 / dp sin^ /3 cos^^^^" /3 / d-f sin^ 7cos3^-i3 7 



47r2 Jo 

/ d9^ sm0^ d9y sinOy I dip f {a 13) g {a ; i = l,2,3. 

Jo JQ Jo 

The angles aij can be determined by p sin aij ~ rjij through the relations 



(C20) 



\/3 1 

m = —m-2 + -j^VN-i , rfii = ^ + ^^'7w-2 + , (C21) 

[211 [21 

mi = Y 2^-3 + ^^^^-2 - 'jP^-^ ' %4 = Y g w-3 - • (^^22) 



3. Results in the 5-limit 

The integrals in the short-range limit, when the range 
h of V{rij) is much smaller than the size scale p, are: 

dr w(a34)(/)(ai3) ~ vi(a)Rl^-^^r>^(j){a) , (C23) 

j dTv{aMaiz) - , (C24) 

j dT v{a3i)c^{a3i) = Ri1'^'^v4>{a) ~ vi{a)4>{Q) , (C25) 

J dr v{a3i)<j){a35) ~ wi(a)^^4^35 0(a) , (C26) 

J dr i;(a34)0(a56) ~ wi(a)J?^^"^V(a) , (C27) 
I 



J dr v(ai3)</)(ai3) = '^^^Z-H =^ «2(a)</)(0) ,(C28) 

('2') 

dr v{ai3)(f>{a^i) ~ U2(a)^i3i4'/'(") ; * = 1,3 , (C29) 

J dr v{ai3)(f>{a23) ^ W2(a)0(a) , (C30) 

/ dr w(ai3)(?!)(a24) — W2(a )i?i|4^(a) , (C31) 

y dr v(ai3)0(a45) ^ i'2(«)i?i345</'(a) • (C32) 



The integrals are given by 



^SL </>(«) = -r rhN-i2\ / rf7sin'7COs3^^''7 , (C33) 

1 ( — 2 — j "'o 



19 



where sinagg = VS cos a sin 7/2, sinQ!4g = cos a cos /3o sin 7, sin/3o = ta,na/\/3, and 



i?g,0(a) EE ,^3^1,,; / d7sin2 7COs3^-i3 7 / d9, sin 9, , 



2 r(M^) r/^ 



2 ) "'O 
2 



sin^ = - sin^ a + - cos^ a cos^ /3o sin^ 7 + sin a cos a cos /3o sin 7 cos ( 



2V2 
3%/3 ' 



4 2 4:\/^ 

sin^ Q!24 = - sin^ a + 7: cos^ a cos^ /3o sin^ 7 H p sin a cos a cos /3o sin 7 cos 

9 3 3v3 



,2 „,0 _ 1 „;„2 „ , 1 „2 ^ „;„2 . 



1 



sm aiQ = — sm a -\ — cos asm 7H ^ sm a cos a sm 7 cos fj- . 

"42 \/2 / ^ 



(C34) 
(C35) 

(C36) 
(C37) 



^ijki't'i'^) appears as a two-dimensional integral but can eq. (B7). 
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